// Replication code for 'Measuring support for European integration using a Bayesian IRT model'
// Author: Michele Scotto di Vettimo - m.scotto-di-vettimo@exeter.ac.uk
// Journal: European Union Politics

cd
use "main_data.dta",clear

set scheme s1mono

// Install required packages
*ssc install grc1leg

// Figure 1: DR algorithm estimates compared to the "membership" question in Germany
foreach v of varlist AHmood GSmood memb_ratio_posneg{
	sum `v' if country_iso2=="DE"
	g z_`v'=(`v'-`r(mean)')/`r(sd)' if country_iso2=="DE"
}
lab var z_memb_ratio_posneg "EU membership: positive/(positive+negative)"
lab var z_AHmood "Anderson & Hecht 'mood'"
lab var z_GSmood "Guinaudeau & Schnatterer 'mood'"
  
line z_* t if country_iso2=="DE" & t>1972 & t<2020, title("(a)",position(11)) ytitle("") xtitle("") xlab(1975 "1975" 1985 "1985" 1995 "1995" 2005 "2005" 2015 "2015") lwidth(medthick medthick medthick) lcol(gs10 gs10 black) lpattern(line dash line) leg(cols(1)) name(gDE1,replace) 

g p1=eumemb_Yes
g p2=eumemb_Yes+eumemb_Eq
g p3=100
lab var p3 "'a bad thing'"
lab var p2 "Neutral/DK"
lab var p1 "'a good thing'"
lab var AHmood "G&S 'mood'"
lab var GSmood "A&H 'mood'"
twoway (area p3 t,color(gs2) lcol(gs2)) (area p2 t,color(gs10) lcol(gs10)) (area p1 t,color(gs14) lcol(gs14) title("(b)",position(11))) || (line GSmood t) (line AHmood t,lpattern(dash)) if country_iso2=="DE" & t>1972 & t<2020,ytitle("Cumulative % 'EU membership is...'") xlab(1975 "1975" 1985 "1985" 1995 "1995" 2005 "2005" 2015 "2015") ylab(0 "0" 20 "20" 40 "40" 60 "60" 80 "80" 100 "100") xtitle("") legend(order(1 4 2 5 3)) name(gDE2,replace)

graph combine gDE1 gDE2, graphregion(fcolor(white)) col(2) name(combined, replace)
graph display combined, xsize(6) ysize(3) 
*gr export "germany.pdf", as(pdf) replace 

// Figure 2: EU support estimated with a Bayesian IRT model. Shaded area indicates ± 2 standard deviations
lab var opinion "IRT model"
levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
twoway (rarea opinion_plussigma opinion_minussigma t if country_name=="`c'" & policy=="jags",  fcolor(gs13) lcolor(gs13)) ///
 (line opinion t if country_name=="`c'" & policy=="jags",title("`c'") ytitle("") xtitle("") xlab(1975 "75" 1985 "85" 1995 "95" 2005 "05" 2015 "15",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick) lcol(black) legend(off) name(g`num',replace))
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames' g0,ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "jags_by_country_with_sigma.pdf", as(pdf) replace 
loc graphnames

// Figure 5: Comparison between the IRT model and other indicators of EU support
twoway scatter bnfit eumemb_Yes unific trustEU euimg t if country_iso2=="EU", mcol(gs10 gs14 gs12 gs12 gs2) msize(small small small small) || line opinion opinion_plussigma opinion_minussigma t if country_iso2=="EU" & policy=="jags", lcol(gs2 gs2 gs2) lpattern(. dash dash) title("") ytitle("") xtitle("") xlab(1975 "1975" 1985 "1985" 1995 "1995" 2005 "2005" 2015 "2015") ylab(20 "20" 40 "40" 60 "60" 80 "80") lwidth(medthick) name(gitems,replace) legend(order(1 "Benefit from EU" 2 "EU membership" 3 "Unification of Europe" 4 "Trust in EU" 5 "EU positive image" 6 "IRT model") cols(3))
*gr export "irt_vs_items.pdf", as(pdf) replace 

// Figure A2: Comparison between IRT (black) and Dyad Ratios (gray) estimates of support for Europe.
levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
corr opinion dr_alg if country_name=="`c'"
loc correlation = substr("`r(rho)'",1,3)
line opinion dr_alg t if country_name=="`c'" & policy=="jags",title("`c'") ytitle("") xtitle("") xlab(1975 "75" 1985 "85" 1995 "95" 2005 "05" 2015 "15",angle(45))  ylab(0 "0" 20 "20" 40 "40" 60 "60" 80 "80" 100 "100") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick medthick) lcol(black) name(g`num',replace) legend(off) text(25 2015 "r = `correlation'")
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames' g0,ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "jags_by_country_vs_dr.pdf", as(pdf) replace
loc graphnames

use policy_areas.dta, clear

preserve
keep if policy=="social"

// Figure 3: Support for EU integration in social affairs by type of welfare model.
line nordic central eastern southern angl balt t, title("") ytitle("") xtitle("") xlab(1990 "90" 1995 "95" 2000 "00" 2005 "05" 2010 "10" 2015 "15" 2020 "20",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") lpattern(. . . dash dash dash) lcol(gs2 gs12 gs8 gs2 gs12 gs8) legend(cols(3)) // Do not plot sigma here as lines overlap too much
gr export "social_by_region.pdf", as(pdf) replace 

// Figure A5: Bayesian IRT estimates of public support for integration in social policies. Shaded area indicates ± 2 standard deviations
levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
twoway (rarea opinion_plussigma opinion_minussigma t if country_name=="`c'",fcolor(gs13) lcolor(gs13)) ///
(line opinion t if country_name=="`c'",title("`c'") ytitle("") xtitle("") xlab(1990 "90" 1995 "95" 2000 "00" 2005 "05" 2010 "10" 2015 "15" 2020 "20",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick) lcol(black) legend(off) name(g`num',replace))
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames' g0,ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "social.pdf", as(pdf) replace 
loc graphnames
restore

preserve
keep if policy=="asylum"

// Figure 4: Support for EU integration of migration policies in six groups of countries.
line nordic central nonvisegrad southern angl visegrad t, ///
title("") ytitle("") xtitle("") xlab(1990 "90" 1995 "95" 2000 "00" 2005 "05" 2010 "10" 2015 "15" 2020 "20",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") ///
lpattern(. . . dash dash dash) lcol(gs2 gs12 gs8 gs2 gs12 gs8) legend(cols(3))
*gr export "asylum_by_region.pdf", as(pdf) replace 

// Figure A6: Bayesian IRT estimates of public support for integration in asylum and immigration policies. Shaded area indicates ± 2 standard deviations.
levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
twoway (rarea opinion_plussigma opinion_minussigma t if country_name=="`c'",fcolor(gs13) lcolor(gs13)) ///
(line opinion t if country_name=="`c'",title("`c'") ytitle("") xtitle("") xlab(1990 "90" 1995 "95" 2000 "00" 2005 "05" 2010 "10" 2015 "15" 2020 "20",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick) lcol(black) legend(off) name(g`num',replace))
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames' g0,ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "asylum.pdf", as(pdf) replace 
loc graphnames
restore

preserve
keep if policy=="econom"

// Figure A7: Bayesian IRT estimates of public support for integration in economic, monetary and fiscal matters. Shaded area indicates ± 2 standard deviations.
levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
twoway (rarea opinion_plussigma opinion_minussigma t if country_name=="`c'",fcolor(gs13) lcolor(gs13)) ///
(line opinion t if country_name=="`c'",title("`c'") ytitle("") xtitle("") xlab(1990 "90" 1995 "95" 2000 "00" 2005 "05" 2010 "10" 2015 "15" 2020 "20",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick) lcol(black) legend(off) name(g`num',replace))
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames' g0,ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "econom.pdf", as(pdf) replace 
loc graphnames
restore

preserve
keep if policy=="market"

// Figure A8: Bayesian IRT estimates of public support for integration in competition, consumer protection and single market rules. Shaded area indicates ± 2 standard deviations.
levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
twoway (rarea opinion_plussigma opinion_minussigma t if country_name=="`c'",fcolor(gs13) lcolor(gs13)) ///
(line opinion t if country_name=="`c'",title("`c'") ytitle("") xtitle("") xlab(1990 "90" 1995 "95" 2000 "00" 2005 "05" 2010 "10" 2015 "15" 2020 "20",angle(45)) ylab(20 "20" 40 "40" 60 "60" 80 "80") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick) lcol(black) legend(off) name(g`num',replace))
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames' g0,ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "market.pdf", as(pdf) replace 
loc graphnames
restore

// Construct validity: predicting eurosceptic vote share
use "election_data.dta",clear

// Models in Table A4: Tobit regression models of public support for Europe and Eurosceptic vote share.

// Model 1
tobit share opinion, ll(0) ul(100) 

// Model 2
tobit share opinion l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, ll(0) ul(100) 

// Model 3
tobit share dr_algorithm l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, ll(0) ul(100) 

// Model 4
tobit share eumembership l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, ll(0) ul(100) 

// Adjust labels to produce coefplot

lab var l_share "Previous vote share"
lab var election "EU election (dummy)"
lab var euid "Exclusive Nat. identity"
lab var fbpop "% foreign population"
lab var fblabfo "% foreign workers"
lab var net_contribution "EU Contribution"
lab var gdp "GDP growth"
lab var une "Unemployment"
lab var t "Year (trend)"

g x=.
replace x=opinion
lab var x "EU support"

tobit share x l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, ll(0) ul(100) 
reg share x l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, r
estimates store m1

replace x=dr_algorithm
tobit share x l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, ll(0) ul(100) 
estimates store m2

replace x=eumembership
tobit share x l_share election euid fbpop fblabfo net_contribution gdp une t i.country_enc, ll(0) ul(100) 
estimates store m3

// Figure 6: Marginal effects of the predictors on Eurosceptic vote share in three models
coefplot (m1,offset(0.2)) m2 (m3,offset(-0.2)), drop(_cons *.country_enc) xline(0, lpattern(dot)) xtitle("Eurosceptic vote share") ///
 p1(label(IRT model) msymbol(O) mcolor(gs2) lcolor(gs2)) ///
 p2(label(DR algorithm) msymbol(D) mcolor(gs8) lcolor(gs8)) ///
 p3(label(Membership) msymbol(S) mcolor(gs12) lcolor(gs12)) ///
 legend(cols(3)) name(cfplot2, replace)
graph display cfplot2, xsize(16) ysize(10)
*gr export "coefplot_complete.pdf", as(pdf) replace 

// ROBUSTNESS CHECKS

// Figure A3: Comparison between main IRT estimates and those calculated excluding trust items

use "main_data.dta",clear
twoway line opinion t if country_iso2=="EU" & policy=="jags" || line opinion t if country_iso2=="EU" & policy=="jags_notrust", xtitle("Year") ytitle("") ylab(30 40 50 60 70) legend(order(1 "All items" 2 "Without trust items"))
*gr export "jags_trust.pdf", as(pdf) replace

// Figure A4: Comparison between IRT estimates using informative prior distributions and those using point priors
use "irt_estimates_point_priors.dta",clear
keep t country_iso2 country_name policy opinion
rename opinion opinion_pointprior
merge 1:1 country_iso2 policy t using "main_data.dta", keep(3)

levelsof country_name,local(countries)
qui foreach c of local countries{
sum country_eu_seq if  country_name=="`c'"
loc num = `r(mean)'
corr opinion opinion_pointprior if country_name=="`c'"
loc correlation = substr("`r(rho)'",1,3)
line opinion opinion_pointprior t if country_name=="`c'",title("`c'") ytitle("") xtitle("") xlab(1975 "75" 1985 "85" 1995 "95" 2005 "05" 2015 "15",angle(45))  ylab(0 "0" 20 "20" 40 "40" 60 "60" 80 "80" 100 "100") yline(50, lpattern(dash) lcol(gs4)) lwidth(medthick medthick) lcol(black) name(g`num',replace) legend(off) text(25 2015 "r = `correlation'")
if `num'!=0{
	local graphnames `graphnames' g`num'
}
}
graph combine `graphnames',ycommon graphregion(fcolor(white)) col(5) name(combined, replace)
graph display combined, xsize(15) ysize(18)
*gr export "jags_compare_priors.pdf", as(pdf) replace